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Abstract 

Pearle (1970) gave an example of a local hidden variables model which exactly reproduced the singlet 
correlations of quantum theory, through the device of data-rejection: particles can fail to be detected in a 
way which depends on the hidden variables carried by the particles and on the measurement settings. If 
the experimenter computes correlations between measurement outcomes of particle pairs for which both 
particles are detected, he is actually looking at a subsample of particle pairs, determined by interaction 
involving both measurement settings and the hidden variables carried in the particles. 

We correct a mistake in Pearle’s formulas (a normalization error) and more importantly show that the 
model is more simple than first appears. We illustrate with visualisations of the model and with a small 
simulation experiment, with code in the statistical programming language R included in the paper. Open 
problems are discussed. 1 


Pearle’s model simplified 

Pearle (1970) considered a source emitting pairs of particles. The particles carry hidden variables X and 
Y which we take to be random points in the unit ball in R 3 . We assume that Y = —X and X/0 with 
probability 1. Write X = RU where ||U|| 2 = 1 and R > 0. One might think of the unit length vector U as 
the direction of spin of the first particle, and —U as the direction of spin of the second, equal and opposite 
points on the unit sphere S 2 , while the scalar R is some kind of amplitude of spin. 

Assume that the direction U is uniformly distributed on S 2 and statistically independent of the amplitude 

R e (o,i]. 

Notation: bold (as opposed to italic) indicates a vector; random vectors and random variables are denoted by 
upper case symbols, while lower case is used for non-random quantities. 

Each particle gets measured in directions a and b respectively (points on S 2 , chosen freely by the experimenter). 
The possible outcomes are +1 (“spin up”), —1 (“spin down”), and last but not least “no detection”, according 
to the following rule: if the angle between X and a is less than Rn/2 then the outcome of measuring the first 
particle is +1; if the angle between X and —a is less than Rn/2 then the outcome of measuring the first 
particle is —1; otherwise the particle is not detected at all. The rule for the second particle is exactly the 
same story as for the first particle, with X and a replaced throughout by Y and b. 

The smaller of the angles between X and ±a is cos -1 |U • a| so the recipe becomes: the outcome of measuring 
the first particle is signU • a if cos -1 |U • a| > Rn/2 while there is “no detection” if cos” 1 |U • a| < Rn/2-, 
the outcome of measuring the second particle is —signU • b if cos -1 |U ■ b| > Rn/2 while it is not detected if 
cos" 1 |U • b| < Rn/2. 

Pearle (1970) gives a formula, (22) in his paper, for a particular choice of the probability density of R (but 
take note of his idiosyncratic normalisation (1)!). There is an error in his derivation, as can be verified 
by integrating the density over the whole range: combining (1) and (22), we get a density which does not 
integrate to 1. Working through Pearle’s paper in detail, it turns out that the only error in (22) is the 
normalisation constant, and this probably derives from an incorrect normalization in (14) where Pearle 
switches from R to S = cos(Rn/2), but it is difficult to be certain about this, since his notion of probability 
density is ambiguous and unconventional. 

1 This paper has been composed in the R markdown language and typeset using the R package "knitr" and RStudio, a popular 

IDE for working with R. The source code contains therefore interleaved passages of text (including LaTeX code, for instance, for 

mathematical formulas) and R code. Processing it in R with knitr generates a LaTeX source file containing interleaved text, R 
code and R textual output. It also generates pdf figures of R’s graphical output. Finally pdflatex generates the pdf you are 
reading now. The source file is available from the author. 
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Here I will present an alternative and much simpler description of the distribution of R and also of the whole 
model, via the distribution of S = cos(f?7r/2) e [0,1). It turns out that the distribution of S can be expressed 
by the formula S = (2/\/V) — 1 where V ~ Unif(l,4); and moreover it is S which we primarily need to know 
in order to simulate the model. 

In terms of S' = ( 2/\/V ) — 1, the recipe for simulating the measurement of one pair of particles is as follows: 
generate U uniformly at random on the sphere S 2 and independently thereof, generate V uniformly at 
random in the real interval [1,4]. Compute S = (2/yfV) — 1 € (0,1), 4 = U-a, and B = — U • b. Particle 1 
is detected if and only if |A| > S, and if it is detected, the outcome of measurement is sign(A). Particle 2 is 
detected if and only if \B\ > S, and if it is detected, the outcome of measurement is sign(H). 

Pearle’s main result is that this model reproduces the singlet correlations: 


E^sign(A)sign(B) |A| > S and 


\B\>S 


— a • b. 


I do not reproduce Pearle’s (magnificent but of necessity very involved) proof. Instead I will just derive 
the density of R according to my specification, so that the reader can compare with Pearle’s formula. I 
will then “prove” Pearle’s result by a simulation experiment. In fact, I would dearly like to see a short-cut 
derivation of Pearle’s result. Through some quite brilliant calculations, he characterizes all possible probability 
distributions of R (equivalently, of S) which will reproduce the singlet correlations as (up to normalization) 
the positive functions within the range of a certain differential operator, and then shows that the operator 
when applied to the constant function the most simple choice one could make - is indeed positive. Further 
details are given in an appendix at the end of this paper. 

According to my definitions, R = cos -1 (S')/ (7t/2) and it follows that for r € (0,1), 

Pr(i? < r) = Pr(S > cos(rf)) = Pr(2/v / P — 1 > cos(rf)) 


= Pr ( W < 


= Pr \ V < 


1 + cos(rf) 
4 


(1 + cos(rf)) 2 
4 


3 l (1 + cos(r|)) 2 


- 1 


and hence the probability density of R is 


f n 4 7T sin(rf) 
jR[r> 3' ‘ 2 ‘ (l + cos(rf)) 3 

47T sin(rf) 


3 (l + cos(rf)) 3 

on the interval (0,1). Compare this to Pearle’s formulas (1) and (22) combined: 


47rp(r)r J = 


2 16 sin(rf) 


3 (l + cos(rf)) 3 ' 


Here is a graph of the probability density of R , together with a rough numerical check that it integrates to 
1. Since the probability density is monotone increasing we get guaranteed lower and upper bounds to the 
integral by summing the value of the density on a regular grid of points between 0 and 1, omitting the right 
hand and left hand endpoints respectively, and dividing by the number of intervals generated by the grid. 
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r <- seq(from = 0, to = 1, length = 1001) 

f <- (4 * pi / 3) * sin(r * pi / 2) / (1 + cos(r *pi / 2))~3 
plot(r, f, type = "1") 



r 


sum (f [1:1000] / 1000) 

## [1] 0.9979072 
sum(f [2:1001] / 1000) 

## [1] 1.002096 

If the points X had been chosen uniformly distributed within the unit ball, the probability density of their 
distance R to the orgin would have had probability density 3r 2 , 0 < r < 1. In the next plot I compare the 
two densities (the one corresponding to a uniform distribution over the ball in green). 

g <- 3 * r~ 2 

plot(r, f, type = "1") 

lines(r, g, col = "green") 
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r 

We see that Pearle’s points have a tendency to be closer to the surface of the ball than if they had been 
uniformly distributed throughout it. 

According to Pearle’s model, a point in the ball is observed when its spin is measured in a certain direction if 
and only if the point lies in one of two, what Pearl describes as “mushroom shaped”, regions around the 
measurement direction and its opposite. In the following plot I draw the intersection of the two mushrooms 
(one blue, the other red) with a plane through the origin containing the measurement direction, which is taken 
to be the direction of the positive z-axis. I superimpose on this plot a sample of 1000 particles distributed in 
a circularly symmetric way about the centre of the unit disk , with distance to the origin distributed as in 
Pearle’s model. The result is a 2D characture of Pearle’s 3D model: some statistical features are the same, 
some are different. 

The picture is neither a 2D section nor a 2D projection of the 3D model. However, it should help the reader 
to visualise the model. The points are coloured blue, red or black according to whether the corresponding 
particle measurement result is an outcome spin up, spin down, or the particle is not detected. 

The actual detection regions (the red and the blue mushroom) are formed by rotating the 2D boundaries 
about the a;-axis. The actual distribution of the 3D hidden variable of the particle being measured has the 
same radial component as in the plot, but its direction is now uniform over the sphere, instead of the circle. 
Thus the 3D density of points is less than what the picture suggests as we move further from the origin, 
however it still increases as we move outward relative to a uniform density. 


set.seed(1234) 
par(pty = "s") 

plot(r * cos(r * pi / 2) , r * sin(r * pi / 2) , type = " 
xlim = c(-l, 1), ylim = c(-l, 1), asp = 1, col = " 

lines(r * cos(r * pi / 2) , - r * sin(r * pi / 2) , col = 

lines(- r * cos(r * pi / 2) , r * sin(r * pi / 2) , col = 

lines(- r * cos(r * pi / 2) , - r * sin(r * pi / 2) , col 

t <- seq(from = 0, to = pi, length = 1000) 
lines(sin(t) , cos(t), col = "blue") 
lines(- sin(t), cos(t), col = "red") 

S <- 2 / sqrt(l + 3 * runif(1000)) - 1 
R <- acos(S) / (pi/2) 


1 ", 

blue" ) 
"blue") 
"red") 

= "red") 
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Theta <- runif(lOOO) * 2 * pi 

Sign <- (Theta < pi/2 I Theta > 3 * pi/2) 

Dev <- Theta * (Theta < pi/2) + (2 * pi - Theta) * (Theta > 3 * pi/2) + 
abs (Theta - pi) * (Theta >= pi/2 & Theta <= 3 * pi/2) 

Obs <- Dev < R * pi/2 
Col <- rep("black", 1000) 

Col [Obs & Sign] <- "blue" 

Col [Obs & ISign] <- "red" 

points (R * cos (Theta), R * sin(Theta) , pch = col = Col) 



r * cos(r * pi/2) 

Pearle found his model by fixing the mushroom shape first, then looking for a probability distribution of the 
radial distance R such that the pair (mushroom shape, distribution of R) reproduce the singlet correlations. 
If we transform the unit ball onto itself in a continuous way by applying a monotone transformation of the 
unit interval onto itself to the distances of points from the origin, we can transform the distribution of R into 
any distribution on [0,1] with cumulative distribution function which is continuous and strictly increasing 
throughout the interval. The mushroom shape will be transformed correspondingly. The author of this paper 
has not discovered a transformation which simultaneously makes both the the mushroom shape and the 
distribution of R more simple than what they are at present. In fact, Pearle’s choice does amount to fixing 
one of these two coupled parameters so that it has a direct physical interpretation: a particle pair with a 
particular value r of R is such that each particle is not detected at all if the direction of its spin, thought of 
now as an undirected line through the origin, deviates by more than r7r/2 from the direction in which the 
spin is measured, also thought of an undirected line through the origin. 

Altogether, Pearle’s derivation of his model was a tour-de-force in imagination, analysis and geometry. 
Whether or not there is a short-cut to getting his results and whether or not they can be improved are 
interesting challenges. As we will see in the next section, the model has one major defect, namely the rate at 
which a pair of particles are both detected depends quite strongly on the pair of settings with which they are 
measured. This phenomenon would be experimentally observable; conversely, the usual quantum mechanical 
modelling of this experiment, and assuming that particles are detected independently of the direction in 
which their spin is measured, predicts that the rate of pair detection is independent of measurement settings. 
So we are left with the open problem: is there a distribution of R reproducing the singlet correlations which 
does not have this defect? Pearle does not answer this question explicitly but his text suggests that he 
believes the answer is negative. A numerical analysis (see Appendix) confirms. 
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A simulation experiment 


We now present a simulation of the model in the statistical programming language “R”. First of all, we (re)set 
the random seed, for reproducibility. To see results based on a fresh sample, replace the (integer) seed by 
your own, or delete this line and let your computer dream up one for you (it uses system time + process ID 
to do this job). 

# set.seed() # Initialise random seed from system time + process ID 
set.seed(9875) # Initialise random seed deterministically 


We will generate uniform random points on sphere generated using the “trig method” (method 3) of Dave 
Seaman: see http://rpubs.com/gillll09/13340 for an R illustration. This very effective but little known 
method uses the coincidence that in 3D, a uniform point on the sphere has a z coordinate which is uniformly 
distributed between —1 and +1. So we proceed as follows. 

(a) Choose Z uniformly distributed in [—1,1]. 

(b) Choose 0 uniformly distributed on [0,27r). 

(c) Let R= y/(l - Z 2 ). 

(d) Let X = Rcos(0). 

(e) Let Y = i?sin(0). 

In the following simulation, the measurement directions will all be in the equatorial plane, so only Z and X 
have been generated and are treated as X and Y. 

First of all we set up the measurement angles for setting “a”: directions in the equatorial plane. 

singles <- seq(from = 0, to = 360, by = 1) * 2 * pi/360 
K <- length(angles) 

corrs <- numeric (K) # Container for correlations 

Ns <- numeric (K) # Container for number of states 

For setting “b” we’ll use a fixed direction. 

beta <-0*2* pi/360 # Measurement direction 'b ' fixed, in equatorial plane 

Then the sample size (number of pairs of particles). 

M <- 10~6 

I use the same, single sample of M = 10 6 realizations of hidden states for all measurement directions. 

z <- runif(M, -1, 1) 
t <- runif (M, 0, 2 * pi) 
r <- sqrt(l - z~2) 
x <- r * cos(t) 

e <- rbind(z, x) # 2 x M matrix 

The M columns of e represent the x and y coordinates of M uniform random points on the sphere S 2 . 
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U <- runif(M) 

s <- (2/sqrt (3*U+1)) - 1 # Pearle's "r" is arc cosine of "s" divided by pi/2 

b <- c(cos(beta), sin(beta)) # Measurement vector 'b' 

Loop through measurement vectors “a” (except last = 360 degrees = first): 

for (i in 1 :(K - 1) ) { 
alpha <- angles [i] 

a <- c(cos(alpha), sin(alpha)) # Measurement vector 'a' 
ca <- colSums(e * a) # Inner products of cols of 'e' with 'a' 
cb <- colSums(e * b) # Inner products of cols of 'e' with 'b' 
good <- abs(ca) > s & abs(cb) > s # Select the 'states' 

N <- sum (good) 

corrs [i] <- sum(sign (ca[good]) * sign(cb [good]))/N 
Ns [i] <- N 

> 

corrs[K] <- corrs[1] 

Ns [K] <- Ns [1] 


Now we are ready to make some plots of the results. 

plot(angles * 180/pi, corrs, type = "1", col = "blue", 
main = "Two correlation functions", 
xlab = "Angle (degrees)", ylab = "Correlation") 
points (angles * 180/pi, corrs, col = "blue", pch = cex = 2) 

lines(angles * 180/pi, cos(angles), col = "black") 

points (angles * 180/pi, cos (angles), col = "black", pch = cex = 2) 

legend(x =0, y = 0, legend = cO'Pearle", "cosine"), text.col = c("blue", 
"black"), lty = 1, col = c("blue", "black")) 

Two correlation functions 



Angle (degrees) 
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In the second plot we zoom in on just part of the curve. 


plot(angles * 180/pi, corrs, type = "1", col = "blue", xlim = c(0, 90), ylim = c(0, 

1), main = "Two correlation functions", xlab = "Angle (degrees)", ylab = "Correlation") 
points (angles * 180/pi, corrs, col = "blue", pch = cex = 2) 

lines(angles * 180/pi, cos(angles), col = "black") 

points (angles * 180/pi, cos (angles), col = "black", pch = cex = 2) 

legend(x = 0, y = 0.4, legend = c("Pearle", "cosine"), text.col = c("blue", 

"black"), lty = 1, col = c("blue", "black")) 

Two correlation functions 
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Angle (degrees) 

Here is an even closer look at part of the curves. 

plot(angles * 180/pi, corrs, type = "1", col = "blue", 
xlim = c(0, 50), ylim = c(0.8, 1), 
main = "Two correlation functions", 
xlab = "Angle (degrees)", ylab = "Correlation") 
points (angles * 180/pi, corrs, type = "b" , col = "blue", 
lines(angles * 180/pi, cos(angles), col = "black") 
points (angles * 180/pi, cos (angles), type = "b" , col = 
legend(x = 0, y = 0.85, legend = c("Pearle", "cosine"), 
"black"), lty = 1, col = c("blue", "black")) 


pch 


2 ) 


"black", pch = cex 

text.col = c("blue", 


2) 









Correlation 


Two correlation functions 



Angle (degrees) 

Here is a plot of the differences between theory and simulation, with an indication of accuracy. 

plot(angles * 180/pi, corrs - cos(angles), type = "1", col = "blue", 
main = "Difference (red: upper bound to +/- 1 standard error)") 
abline(h = 0, col = "black", lwd = 2) 
lines(angles * 180/pi, l/sqrt(Ns), col = "red") 
lines(angles * 180/pi, -l/sqrt(Ns), col = "red") 
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Difference (red: upper bound to +/- 1 standard error) 
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angles * 180/pi 


max(abs(corrs - cos (angles))) 


## [1] 0.002333192 

Finally, a plot of the proportion of observed particle pairs to emitted pairs, as function of the angle between 
the measurement directions. 

plot(angles * 180/pi, Ns / M, type = "1", col = "blue", 

main = "Rate of detected particle pairs", ylim = c(0, 1)) 
abline(h = (2/3)) 
abline(h = (4/3) * (1 - 2/pi)) 
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Rate of detected particle pairs 



angles * 180/pi 

The two horizontal lines are the maximum and minimum detection rates computed by Pearle: 2/3 and 
4/3(1 — 2/7t) = 0.4845 ... of M respectively. Now, if an experimenter is not using pulsed emission of particle 
pairs but they are being emitted in a continuous fashion according to a Poisson process, then the experimenter 
will have no way of knowing that when neither particle is detected, there was still an emission of a particle 
pair. So the loss of 1/3 of all emitted particle pairs will go unnoticed. But the experimenter will be able 
to see that the rate of double detections depends strongly on the difference between the two measurement 
directions - the maximum rate is more than 4/3 times the smallest. Put another way: the rate at which 
particles are detected at one measurement station with no accompanying detection at the other depends on 
the difference between the two measurement directions. Thus Pearle’s model has some very unsatisfactory 
features: assuming a constant emission rate, the experimenter can see that particles are suspiciously being 
rejected in a way which depends on both the settings. It was only in 2008 that Gisin and Gisin came up with 
a new local hidden variable model for the singlet correlations based on data rejection which possesses all the 
symmetries one would require. Moreover, it is amazingly simple. However it seems further from physical 
interpretation than Pearle’s model. 

But Pearle did more than exhibit one concrete local hidden variable model which reproduces the singlet 
correlations: he also characterizes the class of all distributions of R which does the job. This allows us 
in principle to find out if there is a distribution within the class which leads to a model with all required 
symmetries. I believe the answer is negative (and I believe that Pearle knew this too) but I do not have a 
mathematical proof. Numerical evidence (see Appendix) is very strong and inspection of the numerical result 
might help in constructing a proof. 
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simulation reported here was posted to RPubs in March 2014. As far as I know, nobody has investigated 
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Appendix 

Here I reproduce Pearle’s description of the class of all distributions of R such that the singlet correlation is 
recovered from the measurement outcomes of detected particle pairs. 

Let g be a real function on [0,1] satisfying the symmetry requirement g,(x) = /x(( 1 — x 2 ) 1 / 2 ) for all x. One 
could for instance pick /x arbitrarily on [0, l/\/2] and use the symmetry requirement to determine /x on 
(1/V2,1]. It would probably be wise to impose continuity of the derivative of /x at \/V2. Next compute the 
function h through 


h(x) 


d 2 / z 2 
da; 2 y 1 — a; 2 


J {l-z 2 )in(z)dz-J*(l 



If h is nonnegative and integrable, normalize it to a probability density: this should be the probability density 
of S = cos(I?7r/2). The choice /x = constant delivers the particular distribution of S which Pearle further 
investigates and which we have studied here. 

I am not aware of a simple interpretation of the function /x so its role is hard to understand. It is the result 
of applying a certain differential operator, Pearle’s equation (17) 


fi(x) 


(1 — x 2 ) 1 / 2 d.T 


(l-z 2 ) 2 d / , ' 


4a; dx 


g{x) a 


to the function g defined as the probability of detection of a particle pair, expressed as function of s = cos(r7r/2). 

In principle we can therefore see what happens if we specify g = constant, and giving g(x) = Cx{\ — a; 2 ) 1 / 2 ; 
this results in a candidate for the function h and we only have to find out whether or not h can be normalized 
to a probability density (integrable, nonnegative). I was not able to perform this operation analytically. 
However, numerical integration and differentiation delivers us a candidate h which takes both negative and 
positive values. This provides strong evidence that his model does not include a distribution for R (or 
equivalently S) such that the pair detection probability is independent of the settings. The numerical analysis 
did confirm the theoretical analysis for the case /x = constant, so the author does have some faith in its 
results. 

2 1 consider the unpublished 2002 preprint Thompson and Holstein the best and most complete of Caroline Thompson’s 
presentations of her model; her website contains links to some of her earlier papers which she did manage to get published. 
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The investigation is hampered by several misprints in Pearle’s paper: for instance the power in equation (21) 
should be —3 not +3, and throughout, normalisation constants are not to be trusted. The following very 
naive code “computes” the density of S first in the case of constant /i then in the case of constant g. 

rm(list = IsO) # clear workspace 
eps <- 10~{-9> 

Npts <- 10~4 

z <- seq(from = eps, to = 1 - eps, length = Npts) 
x <- z 

KernelO <- outer (z~2, x~2, "/") 

KernelO [lower. tri(KernelO)] <- 1 
KernelO <- sqrt(l - KernelO) 

mu <- 1 

First <- mean(sqrt(l - z~2) * mu) 

Kernel <- KernelO * mu 
Second <- colMeans (Kernel) 

result <- x~2 * (First - Second) / (1 - x~2) 

dens <- diff(diff (result[1:(Npts - 100 + 2)]) * Npts) * Npts 
s <- x[l : (Npts - 100)] 

# "-100" because of numerical instability at right endpoint 
dens <- dens/mean(dens) 

plot(s, dens, type = "1", ylim = c(0, max(dens))) 

# Note vertical offset 0.02 to separate curves: 

lines(s, 0.02 + (1 + s)~(-3) / mean( (1 + s)~(-3) ), col = "green") 



S 

mu <- z * sqrt(l - z~2) 

First <- mean(sqrt(l - z~2) * mu) 
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Kernel <- KernelO * mu 
Second <- colMeans (Kernel) 

result <- x~2 * (First - Second) / (1 - x~2) 

dens <- diff(diff (result[1:(Npts - 100 + 2)]) * Npts) * Npts 

s <- x[l : (Npts - 100)] 

dens <- dens/mean(dens) 

plot(s, dens, type = "1") 

abline(h = 0) 



S 
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